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INTRODUCTION 


Methods  of  differencing  the  Eulerian  form  of  the  hydrodynamic 
equations  are  investigated.  One  wishes  to  find  a  method  for  machine 
calculations  that  could  be  used  in  hydrodynamic  problems  involving 
large  distortions  of  matter  in  two  space  dimensions.  Boundaries  be¬ 
tween  materials  are  to  be  carried  and  moved  through  the  fixed  Eulerian 
mesh.  Given  a  mesh  at  some  time  and  with  proper  boundary  conditions, 
the  problem  will  be  to  carry  the  values  of  the  quantities  stored  in 
each  mesh  point  forward  in  time  explicitly  to  a  small  time  6t  later. 

The  problem  will  first  be  discussed  in  one  space  dimension  and  then 
for  the  case  of  two  dimensions.  This  numerical  method  is  for  compres¬ 
sible  flow  in  the  presence  of  shocks.  The  IBM  Electronic  Data  Processing- 
Machine,  Type  704,  was  used  for  the  numerical  computations. 

The  hydrodynamic  equations  for  mass  flow,  momentum  change,  and  en¬ 
ergy  change  are,  respectively: 


V  .  pV 


(1) 


(2) 
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p  ®  =  -  V  .  (PV) 


(3) 


•where  p  is  the  density,  V  is  the  material  velocity,  p  is  the  pressure, 
and  E  is  the  total  energy  per  unit  mass.  E  is  the  sum  of  the  internal 
energy  per  unit  mass  (£)  and  the  kinetic  energy  per  unit  mass,  i.e., 

E  =  t  +  |  V  •  V. 


> 
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CHAPTER  I 


DIFFERENCING  IN  ONE  DIMENSION 


Values  of  p.,  V.,  and  g.  will  be  considered  as  the  value  of  p,  V, 

3  J  J 

and  g,  respectively,  at  the  center  of  the  jth  cell  of  the  mesh.  Strictly 

speaking,  the  product  of  p .  and  the  volume  of  the  jth  cell  is  the  mass  in 

J 

the  jth  cell.  The  product  of  V.  and  the  product  of  E.  with  this  mass  are 
the  momentum  and  energy,  respectively,  in  the  jth  cell.  For  a  one-dimen¬ 
sional  space  mesh,  a  simple  method  of  differencing  Eq.  0-)  j 

in  the  jth  cell  would  yield  - 


n+1 


n  5t 
“  pj  =  +  BE 


{(0V)>l/2  '  ‘^l/2} 


where  (pV)j+1y2  is  value  selected  for  (pv)  on  the  right  boundary  of 
the  jth  cell  and  (pV)*?'  -jy 2  is  the  value  on  the  left.  By  n  is  meant  that 
the  quantity  is  to  be  evaluated  from  the  mesh  values  at  time  step  n,  and 
n+1  is  the  value  resulting  from  the  solution  of  Eq.  (4)  at  a  time  &t 
later  (one  time  step  later). 


h 


JM-1  JM 


JMfl 


Figure  1 


Equation  (4)  is  in  the  form  which  we  will  call  conservative.  The 

value  of  the  flow  quantity,  (pV),  on  the  right  boundary  of  cell  j  is  the 

same  value  used  on  the  left  side  for  cell  j+1.  To  illustrate  a  result 

obtainable  from  Eq.  (4),  we  take  a  box  (Figure  l)  with  no  flow  out  either 

end,  i.e.,  (pV)1_1/2  »  (Pv)JMfi/2  “  °*  we  find  from  E<1* 

JM  /  \ 

Z  Sx.lp.  -  p.  )=  0  analytically,  i.e.,  mass  is  conserved  exactly. 
j_l  J \  d  J  / 

It  has  been  explicitly  assumed  that  (pV)^+-|y2  in  cell  j  is  equal  (p v)j_i/2 
for  cell  j+1,  i.e.,  the  flow  term  (pV)  is  single  valued  on  a  cell  boundary. 

If  (PV)1  ±j2  and/or  (Pv)JMfl/2  are  not  zero>  1:116:11  ihe  mass  of  sVstem 

will  be  changing  as  specified  by  the  boundary  conditions  on  these  boundary 
values.  It  may  be  noted  that  Eq.  (4)  is  conservative  regardless  of  how 
(pV).  ±j2  and  (pV)j+]y2  are  evaluated  when  (pv) j+1/2  for  cel1  j  is  e<lual 
(pV) j_l/2  for  cell  j+1. 

Assuming  one  starts  at  j  =  1,  with  suitable  boundary  values  for 

n+1 

(pV)i_i/2  and  (pV)^^^,  the  problem  of  evaluating  p^.  in  Eq.  (4)  be¬ 
comes  one  of  evaluating  (pV) ^+-jy2  in  each  cell  j.  (P^)j_i/2  liave 

already  been  calculated  in  the  preceding  cell  by  the  same  method  used  in 


-8- 
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calculating  (pv)^+^g  in  the  present  cell  and  vill  be  used  in  a  conserv¬ 
ative  manner.  Starting  at  j  =  JM,  the  situation  is  simply  reversed. 

Equations  (2)  and  (5)  are  not  in  the  form  for  conservative  dif¬ 
ferencing."  Remembering  that  for  a  scalar  cp,  =  -^2  +  V  •  Vcp,  Eq.  (2) 
in  one  dimension  becomes 


bV  A  bV 
p  3t  +  pV 


bp 

cix 


Adding  V  to  both  sides  and  using  Eq.  (l). 


d(pV)  _  bp  ,v)  bV  d(pV)  ...  d(p  +  PV2) 
“§t - ^  “  vpv;  -  V  —  - - ^ - 


~Sx  ”  vpv;  '  v  1F“  - - 5F 


Equation  (3)  is  changed  to  conservative  form  similar  to  Eq.  (2)  and  in 
differential  form  gives 

d(pE)  _  _  d(pV  +  pVE) 
dt  5x 

Thus,  in  conservative  difference  form,  Eqs.  (l),  (2),  and  (3)  are 
written 
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and 


n+1 


*r- 


n 


J 


5t 

5x 


i 


(pV  +  pVE)3? 


3-1/2 


This  system  of  equations  will  be  solved  after  finding  a  difference  method 
for  the  flow  quantities  on  the  right  side  of  a  cell  in  the  order  presented. 
That  is,  p3?*1  from  Eq.  (4)  will  be  used  to  solve  for  Y^+1  in  Eq.  (5);  p^+1 
and  V?fl  will  then  be  used  to  solve  for  g^+1  J^E^1  =  E?+1  +  ^  (V^+1)  j  in 
Eq..  (6).  The  solution  of  Eqs.  (4),  (5),  and  (6 )  for  all  mesh  points  is 
called  the  n+1  calculation. 

Two  difficulties  in  solving  the  above  system  of  difference  equations 
are  related  to  their  stability  and  their  behavior  in  attempting  to  repre¬ 
sent  shock  discontinuities .  The  method  of  von  Neumann  and  Richtmyer1  to 
solve  both  of  these  difficulties  is  to  introduce  a  fictitious  bulk-viscos¬ 
ity  pressure,  Q,  and  to  replace  the  material  pressure  in  the  above  equations 
by  the  material  pressure  plus  Q.  This  viscosity  spreads  the  shock  region 
and  thus  turns  the  discontinuity  into  a  region  of  continuity  and  gives 
stability  to  the  other  regions  of  flow  because  of  its  second  order  ef¬ 
fects. 

bV 

We  will  consider  viscosity  pressures  of  the  form  Q  -  -  a  Four 

types  of  particular  interest  will  be  characterized  by 


JRV 


dV 


(7) 


* 
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(8) 


.  5x 

°la  -  ci  T  pC 


(9) 

(10) 

0  gives  the  Richtmyer-von  Neumann1  type  viscosity,  a,  ,  gives  the 
x\V 

Landshoff^  type  viscosity  when  C  is  the  sound  speed,  or  gives  the 
effective  viscosity  of  the  "particle- in- cell"  method.  is  a  type 

viscosity  introduced  in  this  paper  with  two-dimensional  machine  calcula¬ 
tions  in  mind.  0  depends  on  the  temperature,  as  does  C,  hut  does  not 
require  talcing  a  square  root  as  does  C.  is  a  constant  used  to  set 
the  magnitude  of  the  viscosity  term.  C^q  is  also  a  constant.  C^q  is 
determined  by  equating  0^  (with  C]L  given)  and  0LO  in  some  region  and 
at  some  time  for  the  problem  at  hand.  If  the  values  of  p,  p ,  and  C  in 
the  desired  region  and  at  the  desired  time  are  given  by  p,  p,  and  C, 
respectively,  then  CL.  =  C  ^£2.  in  the  piston  problem  to  be  presented, 

L0  _  1  2p 

the  values  of  p,  p,  and  C  were  the  theoretically  expected  values  taken 
from  the  reflected  shock  region. 

A  stability  criterion1  for  the  present  system  which  is  not  uncondi¬ 
tionally  unstable  is  not  the  most  pertinent  consideration  for  selecting 
6t.  Experience  has  shown  that  6t  will  have  to  be  roughly  ten  times 
smaller  than  required  by  the  stability  requirement  to  give  good  accuracy. 


0  .  =  Cn  $  p  |V| 

pic  1  2  H 1  1 


and 


°L0  “  CL0P 
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Let  M  (=  be  the  Mach  number.  For  M  >  1,  8t  will  have  to  be  selected 

in  such  a  manner  that  matter  will  require  roughly  ten  time  steps  to  cross 
a  cell.  For  M  <  1,  the  selection  of  St  is  to  be  made  so  that  the  sound 
wave  takes  roughly  ten  time  steps  to  cross  a  cell. 

Second-order  terms  nay  be  added  to  Eqs.  (l)  and  (3)  which  effectively 
spread  the  shock  and  have  the  same  stabilizing  effect  that  the  term  Q  does 
in  Eq.  (2).  We  will  consider  in  one  dimension  a  mass  diffusion  term  in 
Eq.  (l)  and  a  heat  conduction  term  in  Eq.  (3)*  It  should  be  noted  that 
while  the  replacement  of  p  by  the  material  pressure  plus  Q  results  in  a 
dissipative  term^  in  Eq.  (3),  this  term  does  not  have  the  form  of  heat 
conduction. 

The  quantities  on  the  right  boundary  of  a  cell  in  Eqs.  ( ,  (5),  and 
(6)  will  be  designated  by 


(pV^j+l/2  "  (paVa^  "  Co  Dj+l/2^pj+l  “  pj'} 
(p  +  pw)  .+l/2  =  Pa  +  Q  +  (pbV  •  vd 
(PV  +  pVE) .+l/2  =  (p^  +  Q)  Ve  +  (pcVc)  Ea 
“  C2  Hj+l/2^j+l  " 


(11) 

(12) 


(13) 


where  each  separate  appearance  of  p,  V,  p,  and  E  is  distinguished  by  an 
alphabetic  letter  as  a  subscript.  These  alphabetically  subscripted  terms 
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are  <v  V  V  V  V  V  V  V  <V  V  311,1  Ea*  Each  of  these  tems> 

as  discussed  in  the  next  paragraph,  will  be  evaluated  on  the  «j  +  l/2 
boundary  of  the  cell  j.  Dj+1/2  is  the  Portion  of  the  explicitly  added 
mass  diffusion  term  which  is  to  be  considered  a  function  of  p,  V,  and  £; 
E'+l/2  1s  corresPon^nS  function  in  the  heat  conduction  term.  Cq  is 
the  constant  appearing  in  the  explicitly  added  mass  diffusion  term;  C 2  is 
the  corresponding  constant  in  the  heat  conduction  term.  Q  is  the  explic¬ 
itly  added  viscosity  term  as  designated  by  Eqs.  (7)>  (9)>  or  (l^)* 

We  now  consider  various  differencing  types.  Let  f .  represent  any 

tJ 

single  alphabetically  subscripted  quantity  other  than  the  explicitly 
added  diffusion  terms  represented  in  Eqs.  (ll),  (12),  or  (13)  in  the  jth 
cell.  The  types  of  differencing  of  these  terms  considered  here  will  be 


V/2  -  fj  2  — 


Type  II: 


f . 


0+1/2 


=  f .  if  V„  >  0 
0  T 


J  =  0 


if 


=  if 
0+1 


VT  =  0 


VT  <  0 
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r 


Type  III: 


f*,-,  /o<=  0 


>1/2 


-  (6fj +  5V  -  fj-i)/s 


'  (6f0+l  +  3f]  -  fj+2)/8 


Type  IV: 


"  (“j  +  fJ+i  -  fJ-l)/4 


f j+l/ 2  ^  ° 


"  j+l  +  f j  "  f0+2)/4 


if  VT  >  0 

if  vT  =  0 

if  VT  <  0 

if  VT  >  0 

if  VT  =  0 

if  VT  <  0 


Type  I  gives  a  linear  differencing  scheme.  For  Types  II,  III,  and  IV, 
we  define  the  selected  cell  to  be  j  if  VT  >  0  or  to  be  j  +  1  if  <  0. 
VT,  a  test  flow  value,  is  usually  taken  as  =  (Vj  +  V^+^) .  Thus,  j 
is  the  selected  cell  if  the  flow  is  to  the  right  and  j  +  1  if  the  flow 
is  to  the  left  at  j  +  l/2.  Type  II  indicates  that  one  is  to  use  the 
value  in  the  selected  cell  as  the  value  on  the  boundary.  Type  III  is 
the  result  of  a  quadratic  fit  of  f(x),  centered  in  the  selected  cell, 
using  the  values  in  the  selected  cell  and  the  two  contiguous  cells  and 
evaluated  for  x  on  the  j  +  l/2  boundary.  Type  IV  is  an  extrapolation 
from  the  value  in  the  selected  cell  to  the  j  +  l/2  boundary  using  the 
slope  given  by  the  cells  on  both  sides  of  the  selected  cell.  It  may  be 
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noted  that  all  four  types  of  differencing  schemes  can  be  characterized 
by  three  numbers,  say,  and  If  one  writes 


r 


hfj  +  bf  j+i  +  bf  j-i  if  Y  ^  a 

Si +  h  +  b  T 


f0+l/2  \ 


=  0 


if  VT  =  0 


6lf  j+1  *  S2f3  +  bf  j+2  if  V  <  o 


V 


^i +  ^2 +  b 


T 


then  one  has  the  correspondence 


Type 

I  1 

II  1 

III  6 

IV  k 


1 

0 

3 

l 


0 

0 

-1 

-1 


In  the  present  work,  =  (V\  +  V^+^)  is  used  to  evaluate  Va  from  one  of 

the  above  schemes.  The  resultant  value  of  V  is  used  as  the  test  value 

a 

when  differencing  any  of  the  other  alphabetically  subscripted  terms. 
Obviously,  only  Type  I  differencing  of  pa  will  be  considered,  i.e.,  one 
would  not  set  the  pressure  to  zero  on  a  cell  boundary  if  the  velocity 
were  zero  at  that  point. 
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The  preceding  paragraph  merely  outlines  the  prescription  for  carrying 
out  the  various  differencing  schemes.  At  this  point,  it  seems  appropriate 
■^0  give  a  somewhat  more  detailed  discussion  of  the  significance  of  the 
last  three  schemes  from  both  the  mathematical  and  physical  points  of  view. 

Type  II  differencing  gives  a  diffusion  or  second  order  effect.  For 
VT  >  0  on  both  sides  of  a  cell,  the  derivative  of  f  in  cell  j  will  be  rep¬ 
resented  as 


f  -  f 

df  _  j  3-1 

0 


which  would  normally  be  called  the  derivative  of  f  at  j  -  l/2.  An  expan- 

h2 

sion,  f'(x  +  h)  =  f'(x)  +  h  •  f"(x)  +  |r  f’"(x)  +  •••>  yields  f ' ^2 
-  f i  _  f'.'(x)  when  only  the  first  order  correction  term  is  used.  Cor- 


'3  2 

respondingly,  for  V_  <  0  on  both  sides  of  the  cell,  we  find  f  * 


+  i|£  fj(x).  Thus,  for  VT  o  a  constant,  VT  ^  would  be  given  by 


3+1/2 


=*  f*. 


V  ^  =  V 
VT  ^xT  T 
0 


ft 

0 


2  3 


for  Type  II  differencing.  If  Va  were  a  constant  in  Eq..  (ll)  and  p&  were 

treated  by  Type  II  differencing,  then  a  mass  diffusion  term  would  result 

from  this  differencing.  If  V  were  not  constant  and  also  given  Type  II 

a 

differencing,  a  more  complicated  type  of  second  order  velocity  term  vonld 

also  be  present.  Differencing  in  Eq,.  (12)  and  E&  in  Eq..  (lj)  by  Type  II 

3 

gives  second  order  effects  in  the  pic  method. 
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Type  III  differencing  gives  a  fitting  of  the  derivatives,  i.e.,  for 
VT  >  0  on  both  sides. 


ar  6(fj 


-  fi-i>  +  3C V  -  V  -  Wj.i  -  Vs1 

B&x 


"  (6fj-l/2  +  5f j+l/2  ‘  f j-3/2^8 


If  VT  <  0  on  both  sides,  then  f /o  is  used  in  the  derivative  fit. 


j+3/2 


Type  IV  is  similar  except  that  it  uses  extrapolation  of  the  derivative. 
If  Vrj,  >  0  on  the  right  and  <  0  on  the  left  of  cell  j,  then  Types  III 
and  IV  give 


df  fj+l  "  fj-l 


2Sx 


On  the  other  hand,  if  V^,  <  0  on  the  right  and  V^  >  0  on  the  left  of 
cell  j,  then  Type  III  gives 


df 


Tf  -  f  1  Tf  -  f  1 
_  3  XJ+1  Ij-1  1  XJ+2  XJ-2 

~  2  L  26x  J  "  2  j_  ]+Sx 


and  Type  IV  gives 


df 


Tf  -fT  Tf  -fl 

=  ?  3+1  J"1  i  J+j  J"2 

:  [  2Sx  J  L  5§x  _ 
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When  VT  has  different  signs  on  opposite  sides  of  a  cell,  we  have  a  so- 
called  stagnation  region  (low  velocity  region).  The  last  three  equations 
show  how  in  the  jth  cell  would  be  calculated  from  Types  III  and  IV 
differencing  in  a  stagnation  region.  Since  these  latter  calculations 
use  ordinary  linear  type  differencing  for  the  derivatives, or  weighted 
combinations,  there  are  no  diffusion  effects  from  Types  III  and  IV  dif¬ 
ferencing  in  a  stagnation  region. 

The  first  test  problem  selected  for  machine  calculation  was  that  of 
a  steady  state  infinite  shock  moving  to  the  right  and  in  the  leftmost 
cell  at  time  zero.  This  is  called  the  piston  problem.  The  right  side 
of  the  rightmost  cell  (JM)  was  selected  to  be  a  rigid  wall  so  that  the 
reflected  shock  could  be  studied.  The  application  of  the  Hugoniot  rela¬ 
tions  in  this  problem  and  in  the  fractured  diaphragm  problem  to  be  pre- 

k 

sented  later,  are  given  elsewhere  and  will  not  be  discussed  in  detail 

in  this  paper.  These  theoretical  steady  state  solutions  will  be  given 

on  the  graphs  as  straight  lines.  A  polytropic  gas  will  be  assumed.  We 

set  %  -  be  and  p  =  kpO,  where  b  is  the  specific  heat  at  constant  volume; 

k  is  the  gas  constant,  and  9  is  the  temperature .  We  take  b  =  0.06  and 

k  =  0.04  so  that  y(-  1  +  is  j.  In  front  of  the  shock  the  temperature 

and  velocity  are  zero  and  the  density  is  unity.  The  Hugoniot  relations 

then  give  p  =  4,  0  =  8  i,  and  the  shock  velocity  (V  )  equals  4/3  when 
s  s  p  s 

the  material  velocity  behind  the  shock  is  unity.  Thus,  at  the  left 
boundary  and  to  the  left  of  this  point,  the  conditions  are  given  by 

P-1  =  Pl-l/2  “  hr>  V-1  =  Vl-l/2  “  1}  811(1  °-l  =  °l-l/2  =  8  3  *  The 
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condition  for  reflection  at  the  boundary  JM  +  l/2  is  satisfied  by  setting 

VJJMl  *  -VJM>  »»  ’  ®JM>  “d  PjMfl  "  f JM'  »  be  n0ted  that’  r6Sard- 
less  of  the  differencing  scheme  considered  here,  these  conditions  insure 

no  mass  flow,  (pV)J1,I+1/2  =  °>  and  no  enerSy  flow,  (pV  +  pVE)JM+1^2  =  0, 

out  of  the  system  at  the  JM  +  l/2  boundary.  Also,  all  diffusion  terms 

vanish  except  for  Q  at  JM  +  l/2. 

Figures  2  through  29  show  some  of  the  results  of  various  differencing 
methods  for  the  piston  problem.  In  all  these  examples,  p&  =  p^,  pa  = 

=  p  ,  and  V  =  V,  =  V  .  Some  test  problems  were  run  using  p  ^  p,  and 
leaving  Q  out  of  Eq..  (13),  but  all  such  problems  gave  poorer  results  than 
the  ones  shown.  Poorer  results  were  also  obtained  when  the  mass  flow 
terms  were  not  kept  equal.  We  set  (p  V  )  =?  ( p,  V,  )  =  (p  V  )  by  the  above 
equality  of  the  individual  p’s  and  V's.  Whole  terms,  as  contrasted  to 
individual  factors.,  were  differenced  by  the  above-mentioned  types  without 
any  significant  improvements.  An  example  of  this  sort  would  be  to  take 
(paVa)  =  se-*-ec^  TyPes  II,  III,  or  IV  differencing  for  -this 

whole  term. 

When  the  calculation  for  Q  results  in  a  negative  number,  one  gener- 

ally  considers  that  there  is  a  rarefaction  region,  i.e.,  ^  >  0.  Since 

one  does  not  wish  to  spread  a  rarefaction  region,  the  Q  should  be  cut  off 

(Q  =  0)  if  it  is  negative.  The  word  cut  after  the  viscosity  designation 

on  the  graphs  means  that  a  test  on  ^  was  performed  and  Q  was  set  to 
tvtr 

zero  if  ^  >  0.  The  explicit  heat  diffusion  term,  characterized  by 

(Text  continues  on  page  49.) 
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Figure  29 


C£  /  0,  was  cut  off  by  the  same  criterion  on  the  velocity  gradient.  The 
explicit  mass  diffusion  term,  characterized  by  Cq  /  0,  was  cut  off  if  it 
was  negative.  If  the  word  "cut"  does  not  appear  after  the  term  charac¬ 
terized  by  CQ,  C1,  or  Cg,  then  the  corresponding  term  was  not  cut  off. 

The  theoretical  piston  problem  should  have  no  rarefactions,  but  since 
the  machine  calculations  give  oscillations,  the  use  of  cut  off  would 
give  different  results  for  machine  calculations.  The  time  at  which  the 
various  curves  are  plotted  may  be  calculated  (if  desired)  by  knowing  the 
shock  velocity  (4/3)  and  reflected  shock  velocity  (2/3)  iu  combination 
with  the  plotted  position  of  the  appropriate  straight  line  theoretical 
solutions . 

Figures  2  through  6  show  some  results  for  the  case  when  pa,  V&,  P&, 
etc.  are  each  differenced  by  Type  I  differencing,  i.e.,  no  low  order  dif¬ 
fusion  effects  are  introduced  due  to  differencing.  The  density  in  Figure  2 
is  plotted  at  two  times  to  show  that  the  oscillations  do  not  grow  and  though 
the  calculated  shock  lags  the  theoretical  shock,  the  shock  speed  is  cor¬ 
rect.  The  fact  that  the  shock  speed  is  correct  in  all  the  piston  problems 
presented  seems  to  be  the  result  of  the  conservative  differencing  form  and 
not  the  differencing  types.  Throughout,  we  have  taken  the  shock  front  po¬ 
sition  at  which  the  density  has  a  value  halfway  between  value  preceding 
and  the  value  behind  the  shock  front.  In  Figure  2  the  shock  is  late  by 
about  two  cells  because  of  an  initial  delay.  The  lateness  or  earliness 
of  the  shock  front  is  dependent  on  the  type  differencing  used.  Figures  2 
and  3  used  the  Richtmyer-von  Neumann  viscosity  with  the  magnitude  (C^) 
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differing  by  a  factor  of  six.  Figures  3a  and  3b  show  the  shock  before 
and  after  reflection.  E  and  p  of  the  reflected  shock  are  plotted  to 
show  that  the  oscillations  in  p  and  Q  are  out  of  phase  with  each  other. 
This  phase  relation  suggested  the  type  mass  diffusion  coefficient  se¬ 
lected  in  Figure  5.  Figure  5  shows,  in  an  extreme  case,  the  effects  of 
mass  diffusion  on  the  reflected  shock  density.  Here  the  density  is  low 
by  about  15$.  Figure  6  shows  about  the  best  results  that  were  obtained 
for  the  reflected  shock  from  all  Type  I  differencing.  Similar  results 
could  be  obtained  by  using  the  Landshoff /viscosity.  As  expected,  the 
results  in  'the  reflected  shock  region  (stagnation  region)  were  not 
nearly  as  good  when  a  Richtmyer-von  Neumann  or  pic  viscosity  was  used. 

This  is,  of  course,  due  to  their  dependence  on  velocity  rather  than 
temperature.  The  oscillations  in  Figure  2  which  result  from  not  enough 
diffusion  may  also  be  considered  to  be  the  result  of  not  enough  entropy 
change.  Similarly,  Figure  5  would  appear  to  show  too  much  entropy  change. 

Figures  7  through  29  show  the  results  for  the  various  diffusion  ef¬ 
fects  introduced  by  combinations  of  the  four  types  of  differencing  with 
explicit  heat  conduction  and  viscosity.  By  unbounded  in  Figure  11  we 
mean  that  the  oscillations  in  the  reflected  shock  were  greater  than  the 

■  ?Q 

capacity  of  the  machine,  i.e.,  greater  than  about  lCr^  .  The  incident 
shock  is  fairly  good  in  most  of  the  situations  shown.  The  worst  cases 
shown  occur  when  p  is  differenced  as  Types  I  or  II  for  certain  cases 
(Figures  8,  9,  10,  11,  12,  13,  and  21).  Thus,  finding  the  correct 
scheme  for  producing  a  good  reflected  shock  appears  as  the  more 
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difficult  problem.  For  the  reflected  shock  (stagnation)  region,  one  can 
see  from  the  curves  that  one  needs  both  an  explicit  heat  conduction  and 
an  explicit  viscosity  to  produce  the  proper  entropy  change.  We  will  now 
discuss  the  preferred  scheme,  along  with  the  reasons  for  its  selection, 
and  mention  some  possible  modifications.  Thus,  the  curves  will  not  be 
discussed  individually  but  will  be  referred  to  at  various  times  to  il¬ 
lustrate  conclusions. 

We  wish  to  find  the  best  scheme  for  evaluating  the  terns  in  Eqs.  (ll), 
(12),  and  (13) .  We  have  already  seen  that  we  should  take  pa  =  Pb  and  paVa 
=  0  V  =  p  V  .  We  calculate  p  by  Type  I  differencing.  We  are  left  with 
the  problem  of  selecting  p  ,  V  ,  Vd,  Vg,  and  the  three  explicit  dif¬ 
fusion  terms.  We  will  give  V,  and  Eo  Type  II  differencing  for  several 
reasons.  First,  Type  II  takes  the  least  machine  time  of  the  four  types 
considered.  Second,  the  equations  need  diffusion  effects  for  stability. 
Third,  these  diffusion  effects  are  of  the  right  order  of  magnitude  for 
shocks  and  do  not  result  in  excessive  spreading  of  shocks  (see  Figures  7 
through  29).  We  generally  let  VQ  =  V&,  since  V&  will  have  already  been 
calculated,  though  calculations  show  that  its  selection  is  not  very  im¬ 
portant.  We  are  now  left  with  the  choice  of  the  diffusion  terms  and  the 
mass  flow  (p  V  )  as  the  major  problem.  If  a  stagnation  region  or  region 
of  low  velocity  (vhere  the  pic  type  viscosity,  introduced  by  selecting 
as  Type  II,  goes  to  zero)  is  to  be  present  in  a  problem,  then  one 
should  have  an  explicit  viscosity  and  heat  conduction  term.  The  best 
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type  viscosity  is  then,  that  given  by  Landshoff,  as  may  be  seen  from  Fig¬ 
ures  15,  16,  17 ,  18,  23,  and  2k.  If  one  really  worries  about  machine 
time  and  if  one  can  approximate  the  stagnation  region  sound  speed,  the 


viscosity  crTA  may  be  used.  The  best  heat  conduction  coefficient  found 
in  the  present  work  is  given  by  H.+]y2  =  Pacj+1/2  ^see  FiSures  19  through 
29).  Thus,  if  Ojj^  were  used,  the  values  for  would  have  already 

been  calculated.  It  now  seems  that  a  good  value  for  D.+^y2  might  also 
be  given  by  the  product  of  density  and  sound  speed,  though  this  was  not 
tried.  Mass  diffusion  (see  Figures  5 >  8,  9,  10,  and  21)  terms  left  the 
density  too  low  in  the  reflected  shock  in  all  cases  where  the  oscillations 
were  reduced  with  terms  characterized  by  ^  0  and  ^  0.  Therefore, 
the  mass  diffusion  term  is  taken  out  (Cq  =  0).  Figures  2,  3>  k,  6,  11, 

12,  and  13  show  clearly  the  results  of  differencing  pa  as  Type  I.  Thus, 
if  stagnation  regions  are  to  be  present,  p  and  V  should  be  differenced 
as  Type  III  or  IV.  The  above  diffusion  terms  make  the  density  continuous 
in  any  given  material.  This  fact,  plus  the  fact  that  one  wishes  to  move 
the  mass  as  accurately  as  possible  through  the  mesh,  makes  the  author 
prefer  Type  III  differencing  for  p  and  V  .  Type  IV  takes  less  machine 
time  and  appears  nearly  as  good  (see  Figures  28  and  29)  for  this  problem. 
If  there  is  to  be  no  stagnation  region,  then  the  results  given  by  Figure  7 
appear  to  be  the  simplest  and  fastest  method  of  differencing  for  the  pis¬ 
ton  problem.  To  summarize.  Figure  19  gives  -what  we  call  the  best  scheme 
(see  also  Figures  23,  25,  26,  and  27)  if  the  problem  at  hand  is  to  have 
regions  where  the  Mach  number  is  less  than  unity.  The  terms  characterized 


-52- 


by  C  £  0  and  C  p  ^  0  may  be  set  to  zero  in  a  region  of  high  Mach  number 
to  avoid  excessive  spreading  of  shocks.  However,  this  vac  not  done  in 
any  of  the  test  problems  of  this  paper. 

ye  will  now  present  the  results  of  the  study  of  the  fractured  dia¬ 
phragm  problem.  The  same  gas  constants  as  used  for  the  piston  problem 
are  used  for  both  materials  of  the  diaphragm  problem.  The  best  results 
obtained  are  shown  in  Figure  JO.  It  may  be  noted  that  the  "best"  dif¬ 
ferencing  scheme  for  the  piston  problem  was  used.  If  Type  IV  differencing 
is  used  for  p  and  V  ,  then  one  obtains  exactly  the  same  curves.  If  the 
two  explicit  diffusion  terms  are  not  "cut  off,"  then  the  rarefaction  at 
its  left  side  is  four  cells  ahead  of  where  it  should  be  instead  of  the 
two  cells  as  shown.  In  each  case  the  ultimate  speed  was  correct.  Since 
the  rarefaction  wave  and  the  contact  discontinuity  present  the  two  new 
difficulties  in  the  diaphragm  problem  and  since  the  differencing  scheme 
with  "cut  off"  gives  the  best  rarefaction  wave,  the  main  problem  is  the 
contact  discontinuity.  The  initial  values  at  time  zero  and  the  theoret¬ 
ical  values  are  both  represented  as  straight  lines  (Figure  30). 

The  contact  discontinuity  or  interface  is  contained  in  one  cell. 
Physically,  the  velocity  and  pressure  are  continuous  at  the  contact  dis¬ 
continuity.  Thus,  a  desirable  scheme  would  be  to  carry  one  velocity, 
two  temperatures  or  internal  energies,  and  two  densities  for  the  inter¬ 
face  cell.  Schemes  such  as  this  but  which  carry  only  one  temperature 
for  the  interface  cell  were  tried  and  gave  very  poor  results  in 
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Figure  30A 


Sx  =  i ) 


the  temperature  and  density  curves,  because  the  pressure  tends  to  remain 
constant.  Schemes  which  carry  only  one  density  in  the  interface  cell 
were  not  considered,  since  we  did  not  wish  to  allow  one  material  to 
diffuse  into  another.  Call  the  contact  discontinuity  cell  or  interface 
cell  JC.  From  Figure  31>  we  see  that  we  must  carry  the  volume  of  one  of 
the  materials  in  JC.  We  carry  the  volume  of  material  one,  x^,  and  let  x^ 
be  computed  by  xg  =  Sx  -  x^.  The  problem  is  then  to  solve  two  equations 
like  (4)  (one  for  each  material),  one  equation  like  (5)  (with  both  mate¬ 
rials  considered),  and  two  equations  like  (6)  (one  for  each  material)  to 
advance  the  mesh  values  in  JC  to  the  next  time  step.  After  JC  and  the 
rest  of  the  mesh  have  been  advanced  one  time  step  (the  n  +  1  calculation), 
the  value  of  x  is  advanced  one  time  step.  The  equations  used  to  advance 
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the  mesh  values  will  not  he  given  at  this  point  hut  will  he  given  later 
in  more  general  two-dimensional  form  for  cylindrical  coordinates.  The 
conditions  which  are  prescribed  to  solve  Eqs.  (4),  (5),  and  (6)  in  JC 
are:  (pV)  /0  is  zero  for  material  one;  (pv)..  n  is  zero  for  material 


zero  xur  maoex-rtu.  one;  \,pv;^  \j2. 

two;  the  pressure  in  JC  is  the  average  of  material  pressures  one  and  two 

j;  the  total  mass  must  he  accelerated  in  Eq.  (5);  (pV)j+-jy2 

j-l/2  1S 


fc  lP  +  2P 


JC  2 

(the  work  done  hy  the  pressure)  is  zero  for  material  one;  and  (pv) 
zero  for  material  two  ( see  Eqs.  (32)  and  (33)).  The  terms  of  Eqs.  (ll), 
(12),  and  (13)  are  differenced  as  before  except  for  p  near  JC.  When 
considering  material  one,  if  the  differencing  calls  for  a  value  of  den¬ 
sity  in  a  cell  beyond  the  boundary  of  material  one,  the  value  of  the 
density  used  is  that  of  material  one  in  JC.  Similarly,  the  values  for 
material  two  densities  are  extrapolated.  Differencing  p  as  Type  II  at 
only  the  sides  of  cell  JC  gave  poor  results,  is  discontinuous  at 

the  contact  discontinuity  so  that  we  found  a  small  (less  than  2$)  hump 
in  the  velocity  at  cell  JC.  For  this  reason,  we  set  Q  equal  zero  on 
each  side  of  JC  (Figures  30).  The  heat  conduction  terms  used  for  each 
material  in  its  corresponding  energy  equation  were  also  set  to  zero  on 
the  sides  of  JC  which  resulted  in  little  or  no  change  in  the  mesh  values. 

At  first,  an  attempt  was  made  to  change  x^  to  its  new  value  hy  using 
an  interpolated  velocity  at  x^  to  find  its  change  in  one  time  step.  This 
gave  large,  though  hounded,  oscillations  for  all  the  mesh  values  at  JC, 
and  these  oscillations  moved  along  with  the  propagating  shock  and  rare¬ 
faction.  The  reason  for  this  very  poor  result  can  he  seen  as  follows. 
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As  the  value  of  xg  becomes  small,  the  mass  of  material  two,  m^,  in  cell 
JC  also  becomes  small  and  unless  calculated  very  accurately,  the  two 
values  will  combine  to  give  a  poor  value  of  density,  which  then  gives 
poor  pressures,  which  in  turn  lead  to  poor  velocity  and  energy  values. 


To  eliminate  these  difficulties,  the  following  scheme  was  devised.  As 

“2 

x0  and  m0  approach  zero,  we  wish  the  new  value  of  xQ  to  be  such  that  — 
*  d  2 
will  approach  PJC+1*  Also,  as  and  approach  zero,  we  want  —  to 

approach  pJC  1>  When  both  and  xg  have  values  appreciably  different 


from  zero,  we  wish  to  combine  the  above  two  relations  to  obtain  a  mean 

nb  m 

value.  Let  3L  - -  and  X_  =  - :  then  the  desired  relation  for  x_ 

^  Pjc-1  2  PJC+1  1 

at  the  next  time  step  is  given  by 


or  since  Sx  =  x^  +  x^, 

n+1  "  VPJC+1)  +  (5X  -  ^)(ml/pJC-l) 

*1  - - - 


(14) 


All  values  on  the  right  of  Eq..  ( 14)  are  obtained  from  the  mesh  at  the 
end  of  the  n+1  calculation.  Equation  (l4)  insures  that  the  densities 
in  the  interface  cell  will  be  compatible  with  the  densities  in  the  cor¬ 
responding  non- interface  cells  contiguous  to  JC.  Figure  30  indicates 

„  ,  1 

that  this  method  of  calculating  x^  gives  the  desired  result  for  the 
return  to  Eqs.  (4),  (5)>  and  (6)  and  the  start  of  a  new  time  cycle. 
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There  are  two  methods  which  may  be  used  to  determine  whether  or  not 
the  interface  cell  should  be  changed  after  Eq.  (14)  has  been  calculated. 
For  the  material  with  the  smallest  volume,  one  would  test  to  see  if  the 
mass  or  volume  of  this  material  is  likely  to  go  to  zero  in  the  next  time 
cycle.  For  example,  let  x^+1  <  -g-  •  Die  change  in  volume  in  one  time 
step  is  given  by  6x]L  =  x^+1  -  x£.  If  5x]L  >  0,  then  is  not  decreasing 
with  time  and  there  should  be  no  chance  of  JC  changing.  If  5x^  <0,  then 
x^  is  decreasing  and  there  is  a  possibility  of  JC  changing.  Thus,  for 
<  0,  if  (xj+1  +  1.05  6x1)  <  0,  the  cell  should  change  and  if  >  0, 
the  cell  should  not  change.  The  factor  1.05  is  used  only  as  an  example. 
The  excess  over  one  of  this  factor  is  included  to  account  for  accelera¬ 
tions.  To  test  for  a  cell  change  by  using  m^(x^+^  <  ~),  the  most 
appropriate  method  is  to  calculate  by  Eq.  (4)  for  material  one  if  all 
of  mass  one  will  flow  out  of  JC  when  we  go  to  start  a  new  cycle.  Sim- 
ilar  statements  apply  when  x^  >  —  for  material  two.  Die  above  test 
on  mass  was  made  in  the  diaphragm  problem.  The  test  on  volume  takes 
less  machine  time  and  is  more  appropriate  for  two  dimensions.  Equa¬ 
tion  (14)  is  the  same  for  two  dimensions  where  x^  becomes  the  volume 
of  material  one,  Sx  becomes  Idle  total  volume  of  the  interface  cell, 
p  becomes  the  density  of  the  nearest  non- interface  cell  of  type 
two  material,  and  pJC  1  becomes  the  density  of  the  nearest  non- interface 
cell  of  type  one  material.  When  changing  interface  cells  in  one  or  two 
dimensions,  one  applies  the  conservation  of  mass,  momentum,  and  energy 
to  each  material. 
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n+1 

The  calculation  of  x^  represents  a  volume  and  density  change 
since  there  is  no  mass  flow  in  this  boundary  calculation.  Since 
x^  +  x2  =  8x  =  a  const. ,  Sx^  =  -SXg.  Let  p  be  the  pressure  at  the 
interface.  The  work  done  on  material  one  is  -pSx^  and  the  work  done 
on  material  two  is  +pSx1.  Let  and  be  the  internal  energy  per 
unit  mass  in  JC  for  material  one  and  material  two,  respectively.  We 


pbx^  pSx^ 

thus  replace  by  - and  by  +  m^" '  in  bhe  interface  cell. 

In  the  diaphragm  problem  (Figure  30 ),  p  was  calculated  as  the  pressure 
/  1 P  +  \ 

in  JC  (pjC  =  - g - J  using  the  values  of  and  which  were  in  the 

mesh  after  the  completion  of  the  n  +  1  calculation.  Since  the  flow  is 
to  the  right  in  the  diaphragm,  the  values  of  and  are  high  and 
low,  respectively,  at  the  end  of  the  n  +  1  calculation.  The  resultant 
inaccuracy  in  p  gives  the  small  variations  of  the  calculated  density 
and  temperature  near  the  interface  (Figure  30) •  Though  it  has  not 
been  tried,  it  appears  that  the  calculation  of  p  as  p 


-  PJC+1  +  PJC-1 


would  resolve  this  difficulty  in  the  internal  boundary  calculation.  Of 
course,  the  value  of  p  obtained  in  the  n  +  1  calculation  could  be  saved 
and  used  for  the  value  of  p.  This  was  done  in  a  similar  calculation  and 
found  to  give  good  results.  If  p  is  a  given  function  of  space  and  time 

on  an  external  boundary,  then  the  above  method  will  allow  the  correct 

* 

amount  of  work  to  be  done  on  the  system. 

A  difficulty  in  the  diaphragm  problem  (Figure  30)  may  be  noted. 

The  shock  front  is  about  six  cells  wide  as  compared  to  approximately 
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three  cells  for  the  infinite  shock  in  the  piston  problem.  Ihis  indl 
cates  the  selected  difference  scheme  will  have  difficulty  following 
weak  shocks . 

A  cycle  consists  of  the  n  +  1  calculation  and  the  boundary 
calculation.  The  boundary  calculation  consists  of  a  volume  change, 
the  possibility  of  a  cell  change,  and  the  internal,  energy  changes. 
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CHAPTER  II 


DIFFERENCING  IN  TWO  DIMENSIONS 


Equations  (l),  (2),  and  (3)  in  cylindrical  coordinates  with  no 
angular  dependence  may  be  written  as 


dp  _ 

d(rpu)  d(pv) 

(15) 

rdr  dz 

d(pu) 

dp  d(rpuu) 

d(pvu)‘ 

(16) 

“ST“ 

u 

'S 

is 

I 

II 

dz 

d(pv) 

“St 

dp  d(rpuv) 

"Sz  rdr 

d(pw) 

dz 

(17) 

d(pE) 

d(rpu)  d(pv) 

d(rpuE) 

d(pvE) 

(18) 

bt 

rSr  dz 

r^r 

dz 

where  u  and  v  are  the  r-  and  2- components  of  the  velocity. 

We  will  difference  these  equations  for  the  mesh  represented  in 
Figure  32.  The  radius  at  the  center  of  each  cell  is  given  by  r^ 

=  Sr(i  -  l/2)  and  at  the  sides  by  ri+1y2  =  i8r  811(1  ri-l/2  ~  ^  ~  1)5r* 
In  this  notation,  the  volume  and  areas  of  the  cell  (i,j)  are  given  by 
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from  which 


•where  all  mesh  values  on  the  right  are  at  time  step  n.  By  Eq.  (20), 
Eq.  (2l)  may  be  written  as 


"  pi  =  J  [Ai-l/2  pi-l/2  Ui-l/2  "  Ai+l/2  pi+l/2  Ui+l/2 


. j+l/2  j+l/2  j+l/2 

1  pi  1 


or  simply 
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(Apu)i+l/2  +  ^v)J"l/2  -  (Apv)fl/2]  (22) 


Another  method  to  obtain  Eq.  (22)  is  to  integrate  Eq.  (l)  over  the  vol¬ 
ume,  V . ,  convert  the  right  side  by  Gauss'  theorem  to  a  surface  integral, 
and  approximate  the  result  assuming  is  small.  This  latter  method  of 
evaluation  of  the  divergence  is  also  easy  when  other  coordinate  systems 
are  to  be  used.  To  simplify  the  notation,  we  will  designate  quantities 
at  i,  j  -  1/2  as  side  one,  at  i  -  l/2,  j  as  side  two,  at  i,  J  +  l/2  as 
side  three,  and  at  i  +  l/2,  as  side  four.  Thus,  Eq.  (22)  is  written 
as 


.n+1  .n  R.  r  i 

pi  -  pi  =  [(Apu)2  -  (Apu)4  +  (Apv)1  -  (Apv)5J  (25) 

^i 

or  if  mass  is  carried  in  each  cell  instead  of  density 


where  M?  =  p!?  V^.  The  (Apu)  or  (Apv)  terms  are  called  the  mass  flow 
terms.  They  are  the  mass  per  unit  time  flowing  across  their  respective 
sides  of  the  cell.  If  one  starts  at  i  =  1  and  j  =  1  with  suitable  bound¬ 
ary  conditions,  then  we  need  to  calculate  values  for  the  flow  quantities 
only  on  sides  three  and  four.  Other  similarities  to  the  one- dimensional 
cases  are  used,  i.e.,  values  on  sides  one  and  three  are  differenced  as 
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in  one  dimension  in  the  j  direction  for  each  i^  and  values  on  sides  two 
and.  four  are  differenced  in  the  i  direction  for  each  j.  In  Eq.  {2k), 
mass  is  seen  to  he  conserved  as  in  the  one- dimensional  case.  Equa¬ 
tions  (l6),  (17),  and  (18)  are  differenced  similar  to  Eq.  (15)  and 
give 


.n+1  .n+1 

Pi  u. 


.n+1  .n+1 

pJ  v^ 


.n  .n  R+  fV^ 

>i  "i  ’  ^  [si  tP2  ‘  V  +  (APU>2  “2 

-  (Apu)k  u4  +  (flpv^  -  (Apv)5  Ujl 

.n  .n  .  f 

-  pJ  v?  .55 

1  wo  L 


v°. 

5I  (Px  -  Pj)  +  (AP^)2  v2  -  (A^UV  V4 


+  (Apv)1  v1  -  (Apv)3v^J 


(25) 


(26) 


and 


.n+1  .n+1  .n  .n  r 

p?  E1?  -  p?  E?  =  “T  (Apu)  -  (Apu),  +  (Apv)  -  ( Apv ) . 

11  11  -*J  h  ***  - 


x 


(Apu)2  E2  -  (Apu)^  +  (Apv)1  E1  -  (Apv)^  E^J  (27) 


In  the  two-dimensional  example  to  he  presented  below,  the  differencing 
is  like  the  scheme  selected  for  the  diaphragm  problem.  The  p's,  v's, 
and  u's  for  mass  flow  terns  were  calculated  hy  Type  III  differencing, 
p^  is  the  sum  of  the  material  pressure  on  side  four  and  the  "cut  off" 
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value  of  Q  computed  with  the  u's  on  side  four.  Similarly,  is  a  sum 

on  side  three  using  the  v's  in  the  velocity  gradient  portion  of  the 

viscosity  calculation.  The  pressure  times  velocity  terms  in  Eq.  (27) 

are  computed  as  the  sum  of  the  heat  conduction  term  and  of  the  above 

pressures  times  the  velocities  used  in  mass  flow.  This  sum  is  then 

multiplied  by  the  proper  area  so  that  the  equivalent  differential  term 

which  has  been  added  to  Eq.  (3)  is  V  •  CL  D  V£.  The  u's,  v's,  and  E's 

2  E 

which  multiply  the  mass  flow  terms  are  differenced  as  Type  II  with  the 
corresponding  mass  flow  term  as  the  flow  test  number,  V^.  Q  is  zero 
at  r  =  0. 

The  two-dimensional  steady  state  test  problem  geometry  is  shown  in 
Figure  33.  We  wish  to  determine  the  axially  symmetric  flow  about  a  flat¬ 
nosed  cylinder  within  a  cylindrical  pipe.  The  flow  of  the  polytropic 
gas  enters  on  the  left  side  with  a  Mach  number  (m)  of  I.58,  density  and 
sound  speed  equal  1.0,  and  7  =  1.4.  All  cells  to  the  right  of  j  =  14 
were  assumed  to  have  the  same  mesh  values  as  those  with  the  correspond¬ 
ing  value  of  i  in  cell  j  =  14.  The  pipe  and  obstruction  had  rigid  walls. 
Unfortunately,  the  code  was  written  for  the  fixed  configuration  of  Fig¬ 
ure  33  and  for  Sz  =  Sr  =  1.0.  The  left  boundary  was  not  far  enough 
upstream  to  prevent  perturbations  near  the  input  (see  Figures  33,  34, 

35,  36,  and  37)  close  to  the  axis.  Figures  34,  35,  36,  37,  and  38  show 
the  input  conditions  which  were  held  fixed  in  the  j  =  0  cells.  The  flat 
line  in  j  =  0  was  drawn  to  indicate  that  the  input  values  of  p,  u,  v,  and 
0  were  used  for  the  corresponding  value  on  the  right  boundary  of  this 
cell. 
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Considering  these  perturbations  and  the  coarseness  of  the  mesh,  it 

seems  remarkable  that  certain  average  results  were  obtained.  The  posi- 
5 

tion  of  the  detached  shock  is  represented  by  the  solid  line  in  Figure  35* 
The  limits  of  the  force  on  the  flat- nosed  cylindrical  face  are  387*0  and 
33 6.4.  This  gives  a  most  probable  experimental  value  of  361.70.  The  cal¬ 
culated  value  at  steady  state  (Figure  40)  for  this  problem  gave  a  corre¬ 
sponding  value  of  354.0,  i.e.,  about  2$  low.  In  Figures  35  to  39>  the 
mesh  values  at  a  time  of  about  40  are  plotted.  The  steady  state  seemed 
to  be  reached  at  a  time  of  about  20  (see  Figure  40),  which  required  about 
45  minutes  calculation  time.  The  calculation  to  the  time  of  about  40  was 
performed  to  see  if  these  perturbations  of  the  mesh  values  near  the  input 
would  cause  any  change  in  the  steady  state  values.  No  change  was  observed. 
As  an  admittedly  rough  criterion  for  locating  the  shock  position,  we  may 
take  the  position  of  one-half  height  on  the  density  curves  (Figure  34) . 
They  are  plotted  in  Figure  35*  The  perturbation  near  the  input  is  evi¬ 
dent.  The  coarseness  of  the  mesh  and  the  presence  of  diffusion  terms 
seem  to  cause  the  shock  and  reflected  shock  at  the  outside  wall  to  lose 
their  respective  identities  and  to  produce  an  average  position.  In  an 
effort  to  find  some  quantity  which  might  serve  to  differentiate  between 
the  shock  and  reflected  shock,  we  have  plotted  in  Figure  39  "the  tangent 
angle  of  the  velocity  vector.  Since  no  sudden  change  in  slope  is  ob¬ 
served,  it  is  evident  that  this  quantity  does  not  serve  the  purpose. 
Unfortunately,  no  other  single  quantity  examined  has  proved  useful  in 

this  connection.  It  should  be  noted,  however,  that  despite  the  crudity 

(Text  continues  on  page  78.) 
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Figure 


of  the  one-half  height  criterion,  the  position  of  the  shock  front  is  cor¬ 
rect  to  within  a  cell  size.  The  calculated  sonic  line  in  Figure  33  was 
obtained  from  plots  of  the  Mach  number  (Figure  37)*  The  vectors  in  Fig¬ 
ure  33  show  the  direction  of  flow  but  not  the  magnitude  (see  Figures  36 
and  38).  The  normal  component  of  the  velocity  in  the  cells  next  to  the 
flat- nosed  surface  extrapolates  to  zero  (Figure  36).  The  spread  of  about 
six  cells  in  this  weak  shock  is  again  evident  as  it  was  in  the  diaphragm 
problem. 

For  a  two-dimensional  problem  which  has  a  moving  interface  between 
materials  (Figure  4l),  the  only  new  difficulty  introduced  is  that  of  cal¬ 
culating  the  partial  areas  for  the  mass  flow  and  work  (pV)  terms.  The 
criterion  that  each  interface  cell  must  have  two,  and  only  two,  adjoining 
interface  cells  leads  to  six  types  of  three  interface  cell  combinations 
(Figure  4l).  To  maintain  the  interface  normal  to  outside  boundary  walls 
will  require  special  interface  cell  types.  These  special  cells  should  be 
maintained  such  that  only  one  interface  cell  will  be  next  to  an  outside 
boundary  where  the  interface  terminates.  Two  interface  cells  at  i  =  1 
(with  consecutive  j  values),  for  example,  would  violate  the  above  crite¬ 
rion.  The  conservation  of  mass  for  each  material  will  be  given  by  (see 
Eq.  (24)) 

l.n+1  l.n  [  1  1  1  1~| 

-  M?  =  St  I  (Apu)2  -  (Apu)^  +  (Apv)x  -  (Apv)^  I  (28) 


and 


-78- 


•l 


4 


-79- 


Figure  41 


(29) 


2  .n+1  2.n  2  2  2  2 

Mj?  -  t-rf  =  St  (Apu)2  -  (Apu)^  +  (Apv)1  -  (Apv)^ 


vhere  the  superscripts  1  or  2  indicate  the  material.  As  an  illustration 

of  the  evaluation  of  the  mass  flow  terns  in  Eqs.  (28)  and  (29),  consider 

a  Type  I  interface  cell  in  (i,  j)  with  material  1  to  the  right  and  above 

and  with  material  2  to  the  left  and  below.  In  this  case,  A-j^  =  0,  Ag  =  0, 
2  2 

=  A1,  and  Ag  =  Ag.  A  scheme  must  now  be  devised  to  calculate  one  of 

the  partial  areas  on  side  three  and  side  four  since  A,  =  A,  +  A,  and 
12  ?  p  y 
A^  =  A^  +  A^.  The  simplest  scheme  appears  to  be  given  by  using  a  linear 

combination  of  the  partial  volumes,  i.e.. 


2  2 . 

v?  +  v?+1 

1  1 


V<r  -1  - 


-  vrvv 


vi  +  vi+i 


Other  schemes  which  use  the  given  partial  volumes  in  several  cells  and 
attempt  curve  fits  lead  to  quadratic,  cubic,  and  higher  order  equations 
■which  should  be  avoided  if  possible.  A  code  is  being  prepared  which  uses 
the  simple  scheme  mentioned  above.  This  code  will  calculate  the  values 
of  the  density  on  the  cell  boundary  which  joins  two  interface  cells  by 
Type  I  differencing  and  will  use  Type  III  differencing  with  extrapolation 
elsewhere.  The  values  of  u  and  v  are  continuous  through  the  interface 
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and  are  to  be  calculated  as  before.  Pressures  are  calculated  as  in  the 
diaphragm  problem.  The  viscosity  and  heat  conduction  terms  are  to  be 
zero  at  the  boundary  of  any  interface  cell.  Equations  (25)  and  (26)  for 
an  interface  cell  become 


l.ml  2.ml\  .n+1 

M?  +  mJ  /  u” 


{(Apu)g  +  (Apu)2^-u2  -  ^(Apu)j^  +  (Apu)^u^ 

^(Apv^  +  (Apv)^!^  -  ^(Apv)^  +  (Apv^u^j 


and 


/l.ml  2  .ml\  n+1  /l.n  2.n\  n  fvj 

\4  +  4  )  4  -  (4  -  4  h i  =  5t[_sl  (pi  - 

+  ^(Apu)g  +  (Apu)2 j-v2  -  |(Apu)^  +  (Apu)^v^ 

+  -^(Apv)1  +  (Apv).^^  -  |(Apv)5  +  (Apv^v^j 


J 


p3} 


(30) 


(31) 


J 


Equation  (27)  must  be  solved  for  each  material;  thus,  it  takes  the  fol¬ 
lowing  form: 


l.n+1  l.m-1  l.n  l.n  f  1  1  1  1  1  1  1  1 

E"  -  M“  E^  «  St  (Apu)2  Eg  -  (Apu)4  E^  +  (Apv)1  E1  -  (Apv)^  E^ 

+  ^(Apu)2  -  (Apu)^  +  (Apv)1  -  (Apv)5j- J  (52) 
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2.n2.n  [*  2  2  2  2  2  2  2  2 

M±  Ei  =  5t  ^(Apu)g  E2  -  (Apu)^  E^  +  (Apv)1  E1  -  (A pv)^  E^ 

f  2  2  2  2  \  "I 

+  ^(Apu)2  -  (Apu)^  +  (Apv)x  -  (Apv)^  H  (33) 

The  cycling  of  this  problem  would  be  as  follows:  One  first  calculates 
the  values  at  time  n  +  1  for  the  whole  mesh  using  Eqs.  (23),  (25),  (2 6), 
and  (27)  for  non- interface  cells  and  Eqs.  (28),  (29),  (30),  (31),  and 
(33)  for  interface  cells.  Then  the  calculations  are  made  for  each  inter¬ 
face  cell  of  the  moving  boundary.  These  calculations  give  the  new  values 
of  the  partial  volumes,  the  interface  cell  changes,  and  the  new  values  of 
the  internal  energy  resulting  from  the  volume  change.  The  cycle  is  then 
repeated.  This  procedure  is  merely  the  two-dimensional  analogue  of  the 
method  used  in  the  diaphragm  problem. 


2.n+l  2.n*l 

M.  E'? 

1  1 
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